!------------------------------------------------------------------------------
! TJU/Department of Mechanics, Fluid Mechanics, Code START
!------------------------------------------------------------------------------
!> @Author: Jianxin Liu shookware@tju.edu.cn
!> @Date: 2024-08-13 21:38:51
!> @LastEditTime: 2024-08-13 21:38:52
!> @LastEditors: Jianxin Liu shookware@tju.edu.cn
!> @Description:
!> @FilePath: /START_NS_3D/src/splitter/module_splitter_lax_friedrichs.f90
!> @Copyright (c) 2024 by Jianxin Liu email: shookware@tju.edu.cn, All Rights Reserved.
!------------------------------------------------------------------------------
submodule(mod_splitter) submod_splitter_lax_friedrichs
   implicit none

contains

   module subroutine SplitLocalLaxFriedrichs(this, q_in, q_l, q_r, chi_x, chi_y, chi_z, chi_face)
      class(type_splitter_local_lax_friedrichs), intent(inout) :: this
      class(type_vector_base), intent(in) :: q_in(:)
      class(type_vector_base), intent(inout) :: q_l(:), q_r(:)
      class(type_vector_base), intent(in) :: chi_x, chi_y, chi_z, chi_face

      class(type_self_vector), pointer :: q_in_ptr(:), q_l_ptr(:), q_r_ptr(:)
      class(type_self_vector), pointer :: chi_x_ptr, chi_y_ptr, chi_z_ptr, chi_face_ptr

      call vector_2_selfvector(q_in, q_in_ptr)
      call vector_2_selfvector(q_l, q_l_ptr)
      call vector_2_selfvector(q_r, q_r_ptr)
      call vector_2_selfvector(chi_x, chi_x_ptr)
      call vector_2_selfvector(chi_y, chi_y_ptr)
      call vector_2_selfvector(chi_z, chi_z_ptr)
      call vector_2_selfvector(chi_face, chi_face_ptr)

      call SplitLocalLaxFriedrichs2(this, q_in_ptr, q_l_ptr, q_r_ptr, chi_x_ptr, chi_y_ptr, chi_z_ptr, chi_face_ptr)

   end subroutine SplitLocalLaxFriedrichs

   module subroutine SplitLocalLaxFriedrichs2(this, q_in, q_l, q_r, chi_x, chi_y, chi_z, chi_face)
      class(type_splitter_local_lax_friedrichs), intent(inout) :: this
      type(type_self_vector), intent(in) :: q_in(:)
      type(type_self_vector), intent(inout) :: q_l(:), q_r(:)
      type(type_self_vector), intent(in) :: chi_x, chi_y, chi_z, chi_face

      type(type_array) :: v_in(size(q_in)), v_l(size(q_in)), v_r(size(q_in))
      type(type_array) :: k_x, k_y, k_z, k_face
      real(dp) :: vv_in(size(q_in)), vv_l(size(q_in)), vv_r(size(q_in))

      integer :: l, lsta, lend, i

      do l = 1, size(q_in)
         call q_in(l)%get_array_read(v_in(l)%p)
         call q_l(l)%get_array(v_l(l)%p)
         call q_r(l)%get_array(v_r(l)%p)
      end do

      call chi_x%get_array_read(k_x%p)
      call chi_y%get_array_read(k_y%p)
      call chi_z%get_array_read(k_z%p)
      call chi_face%get_array_read(k_face%p)

      do concurrent(l=lbound(v_in(1)%p, dim=1):ubound(v_in(1)%p, dim=1))
         do concurrent(i=1:size(q_in))
            vv_in(i) = v_in(i)%p(l)
         end do
         call SplitPoint(this, vv_in, vv_l, vv_r, k_x%p(l), k_y%p(l), k_z%p(l), k_face%p(l))
         do concurrent(i=1:size(q_in))
            v_l(i)%p(l) = vv_l(i)
            v_r(i)%p(l) = vv_r(i)
         end do
      end do

      call chi_x%restore_array_read(k_x%p)
      call chi_y%restore_array_read(k_y%p)
      call chi_z%restore_array_read(k_z%p)
      call chi_face%restore_array_read(k_face%p)
      do l = 1, size(q_in)
         call q_in(l)%restore_array_read(v_in(l)%p)
         call q_l(l)%restore_array(v_l(l)%p)
         call q_r(l)%restore_array(v_r(l)%p)
      end do

   end subroutine SplitLocalLaxFriedrichs2

   pure subroutine SplitPoint(this, v_in, v_l, v_r, k_x, k_y, k_z, k_face)
      implicit none
      class(type_splitter_local_lax_friedrichs), intent(in) :: this
      real(dp), intent(in) :: v_in(:)
      real(dp), intent(out) :: v_l(:), v_r(:)
      real(dp), intent(in) :: k_x, k_y, k_z, k_face

      real(dp) :: rho, u, v, w, T
      real(dp) :: ubar, c
      real(dp) :: pressure, energy
      real(dp) :: q(5), lambda, E(5)

      real(dp), parameter :: eps = 1.0d-6

      rho = v_in(1)
      u = v_in(2); v = v_in(3); w = v_in(4)
      T = v_in(5)
      pressure = rho*T/this%eos%ref%gamma*(this%eos%ref%flg_uref)**2

      ubar = k_x*u + k_y*v + k_z*w
      c = sqrt(T)*this%eos%ref%flg_uref

      lambda = abs(ubar) + c
      energy = T/this%eos%ref%gamma*(this%eos%ref%flg_uref)**2/(this%eos%ref%gamma - 1.0d0)
      q = rho*[1.0d0, u, v, w, energy + 0.5d0*(u**2 + v**2 + w**2)]
      E = [rho, rho*u, rho*v, rho*w, rho*energy + pressure]*ubar

      v_l = 0.5d0*(E + lambda*q)
      v_r = 0.5d0*(E - lambda*q)

   end subroutine SplitPoint

end submodule submod_splitter_lax_friedrichs
